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We  demonstrate  that  the  requirement  of  Galilean  invariance  determines  the  choice  of  H  function  for  a  wide 
class  of  entropic  lattice-Boltzmann  models  for  the  incompressible  Navier-Stokes  equations.  The  required  H 
function  has  the  form  of  the  Burg  entropy  for  D- 2,  and  of  a  Tsallis  entropy  with  q=  1  -(2 /D)  for  D> 2, 
where  D  is  the  number  of  spatial  dimensions.  We  use  this  observation  to  construct  a  fully  explicit,  uncondi¬ 
tionally  stable,  Galilean-invariant,  lattice-Boltzmann  model  for  the  incompressible  Navier-Stokes  equations, 
for  which  attainable  Reynolds  number  is  limited  only  by  grid  resolution. 
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L  INTRODUCTION 

Lattice-Boltzmann  models  of  fluids  [1,2]  evolve  a  single¬ 
particle  distribution  function  in  discrete  time  steps  on  a  regu¬ 
lar  spatial  lattice,  with  a  discrete  velocity  space  comprising 
the  lattice  vectors  themselves.  The  single-particle  distribu¬ 
tion  corresponding  to  lattice  vector  c,  at  lattice  position  x  and 
time  step  t  is  denoted  by  Nfat).  The  simplest  variety  of 
lattice-Boltzmann  models  employ  a  Bhatnagar-Gross-Krook 
(BGK)  operator  [3],  so  that  their  evolution  equation  is 

V,(x+ c,  ,/ + A0  =  Nfat)  +  N -q(x,r )  -  W,(x,r)], 

for  ,b.  Here  b  is  the  coordination  number  of  the 

lattice,  N^(x,t)  is  a  specified  equilibrium  distribution  func¬ 
tion  that  depends  only  on  the  values  of  the  conserved  quan¬ 
tities  at  a  site,  and  r  is  a  characteristic  collisional  relaxation 
time.  Using  the  Chapman-Enskog  analysis,  it  is  possible  to 
show  that  the  mass  and  momentum  moments  of  the  distribu¬ 
tion  function  will  obey  the  Navier-Stokes  equations  for  cer¬ 
tain  choices  of  equilibrium  distribution  [1]. 

The  viscosity  appearing  in  the  Navier-Stokes  equations 
obtained  from  these  models  is  proportional  to  r— j.  To 
lower  this  and  thereby  increase  Reynolds  number,  practitio¬ 
ners  often  over-relax  the  collision  operator  by  using  values 
of  r  in  the  range  (},!].  For  sufficiently  small  r,  however, 
the  method  loses  numerical  stability,  and  this  consideration 
limits  the  lowest  Reynolds  numbers  attainable. 

In  an  effort  to  understand  these  instabilities,  there  has 
been  much  recent  interest  in  entropic  lattice -Boltzmann  mod¬ 


els  [4-6].  These  models  are  motivated  by  the  fact  that  the 
loss  of  stability  is  due  to  the  absence  of  an  H  theorem  [6]. 
Numerical  instabilities  evolve  in  ways  that  would  be  pre¬ 
cluded  by  the  existence  of  a  Lyapunov  function.  The  idea 
behind  entropic  lattice-Boltzmann  models  is  to  specify  an  H 
function,  rather  than  just  the  form  of  the  equilibrium.  Of 
course,  the  equilibrium  distribution  will  be  that  which  ex- 
tremizes  the  H  function.  The  evolution  will  be  required  never 
to  decrease  H ,  yielding  a  rigorous  discrete- time  H  theorem; 
this  is  to  be  distinguished  from  other  discrete  models  of  fluid 
dynamics  for  which  an  H  theorem  may  be  demonstrated  only 
in  the  limit  of  vanishing  time  step  [7]. 

To  ensure  that  collisions  never  decrease  H>  the  collision 
time  r  is  made  a  function  of  the  incoming  state  by  solving 
for  the  smallest  value  rmin  which  does  not  increase  H,  The 
value  then  used  is  T=rmin//c,  where  0</c<l.  It  has  been 
shown  that  the  expression  for  the  viscosity  obtained  by  the 
Chapman-Enskog  analysis  will  approach  zero  as  k  ap¬ 
proaches  unity  [4-6].  Thus,  the  entropic  lattice-Boltzmann 
methodology  allows  for  arbitrarily  low  viscosity  together 
with  a  rigorous  discrete-time  H  theorem,  and  thus  absolute 
stability.  The  upper  limit  to  the  Reynolds  numbers  attainable 
by  the  model  is  therefore  determined  by  loss  of  resolution  of 
the  smallest  eddies,  rather  than  by  loss  of  stability  [6,8,9]. 

In  a  recent  review  of  the  subject,  Succi,  Karlin,  and  Chen 
[10]  have  pointed  out  that  entropic  lattice-Boltzmann  models 
have  three  important  desiderata:  Galilean  invariance,  non¬ 
negativity  of  the  distribution  function,  and  ease  of  determin¬ 
ing  the  local  equilibrium  distribution  at  each  site  and  at  each 
time  step. 

In  this  paper,  we  shall  construct  entropic  lattice- 
Boltzmann  models  for  the  incompressible  Navier-Stokes 
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equations  which  are  Galilean  invariant  to  second  order  in  the 
Mach  number  expansion  of  the  distribution  function  (quasip¬ 
erfect  in  the  terminology  of  Ref,  [10]),  We  shall  show  that 
the  requirement  of  Galilean  invariance  makes  the  choice  of 
H  function  unique.  We  shall  show  that  the  required  function 
has  the  form  of  the  Burg  entropy  [11]  in  two  dimensions,  and 
the  Tsallis  entropy  in  higher  dimensions.  While  the  analo¬ 
gous  problem  for  the  compressible  Navier-Stokes  equations 
is  difficult  and  remains  outstanding,  the  purpose  of  this  paper 
is  to  point  out  that  the  incompressible  case  is  nontrivial  and 
interesting  in  its  own  right. 

Finally,  a  point  of  clarification:  Throughout  this  paper, 
when  we  describe  the  lattice-Boltzmann  model  as  “incom¬ 
pressible,”  we  really  mean  that  it  is  faithful  to  the  Navier- 
Stokes  equations  only  in  the  asymptotic  limit  of  incompress¬ 
ibility.  This  means  that  the  Mach  number  must  scale  with  the 
Knudsen  number,  and  the  fluctuation  of  density  about  its 
mean  must  scale  with  the  Knudsen  number  squared.  Indeed, 
this  is  the  same  sense  in  which  any  quasicompressible  fluid 
model  may  be  said  to  simulate  incompressible  fluid  equa¬ 
tions.  In  this  asymptotic  limit  the  pressure  is  determined  by 
an  elliptic  equation,  and  the  equation  of  state  becomes  irrel- 


II.  EQUILIBRIUM  DISTRIBUTION 

We  consider  a  Bravais  lattice  of  coordination  number  b  in 
D  dimensions.  We  denote  the  lattice  vectors  by  c, ,  where  i 
=  1, . . .  yb ,  and  their  magnitudes  by  c  =  |c/|.  The  restriction 
to  a  single-speed  model  on  a  Bravais  lattice  is  done  solely  for 
the  sake  of  simplicity  in  presentation.  A  future  publication 
will  generalize  the  results  of  this  paper  to  multispeed  lattice- 
Boltzmann  models  [12]. 

We  demand  that  the  lattice  symmetry  group  be  suffi¬ 
ciently  large  that  the  only  fourth-rank  tensors  that  are  invari¬ 
ant  under  its  group  action  are  isotropic.  The  mass  and  mo¬ 
mentum  densities  are  given  by 


m^i 

i  - 1 


that  H  is  altered  only  by  collisions  and  not  by  propagation, 
we  assume  that  the  H  function  is  of  trace  form 


W  =  S  h(Nt), 

i=  1 

where  h'(x)^0  for  x>0.  If  lim,^ '(•*)  = 00  ■  then  the  nor¬ 
mal  derivative  of  H  goes  to  negative  infinity  on  the  polytope 
boundary,  enforcing  the  non-negativity  constraint  [Eq.  (3)]. 
The  purpose  of  this  paper  is  to  demonstrate  that  the  require¬ 
ment  of  Galilean  invariance  uniquely  determines  the  choice 
of  function  h(x). 

In  passing,  we  note  that  our  choice  of  the  form  of  H 
differs  from  that  of  Karlin,  Ferrante,  and  Ottinger  [5],  which 
is  of  the  form  H  =  2-Niln(Ni/Wi),  where  W,  are  speed- 
dependent  weights  (equal  to  the  global  equilibrium  at  zero 
flow).  That  is,  prior  work  has  allowed  weighted  contributions 
to  H  and  found  solutions  for  which  h  has  the  form  of  a 
(relative)  Boltzmann  entropy,  while  the  present  work  as¬ 
sumes  uniform  contributions  to  H  and  finds  solutions  for 
which  h  is  not  a  Boltzmann  entropy.  Both  approaches  are 
capable  of  yielding  Galilean-invariant  hydrodynamics.  A 
more  general  form  for  H  which  will  subsume  both  ap¬ 
proaches  as  special  cases  remains  an  interesting  theoretical 
challenge. 

The  equilibrium  distribution  function  may  be  found  by 
extremizing  H  with  respect  to  N, ,  subject  to  the  constraints 
[Eqs.  (1)  and  (2)] 

d  I  fi  P  \ 

0=—  \H--p - pu  , 

dNj  \  m  m  / 

where  fi/m  and  film  are  Lagrange  multipliers.  We  quickly 
find 


and  so 


0  =  h'(Nj)  —  /x  —  p-Cj, 


Nf>=^(/u  +  /8-c1), 


and 

b 

pu=2  mCjNi,  (2) 

i  =  1 

where  m  is  the  particle  mass  and  u  is  the  hydrodynamic 
velocity  D  vector.  These  D+  1  quantities  must  be  conserved 
in  collisions. 

If  we  regard  ,  for  i=l,...,fc,  as  coordinates  in  a 
^-dimensional  space,  the  conservation  laws  (1)  and  (2)  re¬ 
strict  the  collision  outcomes  to  a  [b  —  (D+  l)]-dimensional 
subspace.  Since  the  conserved  quantities  are  linear  functions 
of  NiX  the  non-negativity  requirement 

0  (3) 

is  satisfied  within  a  compact  polytope  whose  faces  are  given 
by  the  b  equations,  0  for  /=!,...  tb.  In  order  to  ensure 


where  the  function  <f>  is  the  inverse  function  of  h ' .  The  con¬ 
stants  fi  and  p  are  determined  by  Eqs.  (1)  and  (2).  It  is 
usually  difficult  to  find  an  exact  analytic  expression  for  them 
in  terms  of  the  conserved  quantities  p  and  pu,  though  some 
equilibria  are  known  for  which  this  is  possible  [8,9].  Alter¬ 
natively,  one  may  solve  for  them  numerically  or  perform  a 
Taylor  expansion  in  Mach  number.  We  adopt  the  latter  ap¬ 
proach  below. 

III.  GALILEAN  INVARIANCE 

We  seek  to  Taylor  expand  the  equilibrium  distribution  in 
Mach  number  because  (i)  we  can  do  so  analytically,  (ii)  only 
the  first  two  terms  of  that  expansion  determine  the  form  of 
the  incompressible  Navier-Stokes  equations,  and  (iii)  that  ex¬ 
pansion  is  a  useful  initial  guess  for  any  numerical  solution. 
From  general  symmetry  arguments,  it  is  clear  that  P  will  be 
proportional  to  the  hydrodynamic  velocity  u,  so  that  we  may 
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begin  our  Mach  number  expansion  by  expanding  Eq.  (4)  for 
small  p.  We  get 

Ap=*(/i)  +  *'(/*)0-«»+  \<t>"(n)W  c,c,+ 

Inserting  this  into  Eqs.  (1)  and  (2),  and  using  general  prop¬ 
erties  of  the  Bravais  lattice,  we  find 

mbc 2  . 

p=mb<i>(fi)+  (fi)P2+  ■  ■  ■ 


D  U (x)*”(x) 

D  +  2I  [ <t>'(x)]2  ' 

If  the  factor  g  is  not  unity,  Galilean  invariance  will  be  bro¬ 
ken.  Thus  we  demand  g  =  1 ,  and  this  yields  the  second-order 
nonlinear  differential  equation 

^(x)rW=(i  +  |)[<^'U)]2. 

The  general  solution  to  this  equation  is  of  the  form 


and 


<f>(x)  =  CDr2(x —aC)r, 


mbc2 


pu= 


D 


where  the  ellipses  denote  third  or  higher  order  terms  in  Mach 
number.  Inverting  this  perturbatively  we  find  that,  to  second 
order  in  Mach  number,  the  Lagrange  multipliers  are  given  by 


D 


fi=x- 


P_ 

mb 


*"W 


2c2  [<f>'(x)]2 


-«2+. 


where  x=hf(plmb)y  and  by 


P 


0= 


D  mb 

- u_l_ . . . 

c2  <f>'(x) 


Inserting  these  into  Eq.  (4),  we  obtain  the  equilibrium  distri¬ 
bution 


1 


D  D 2  <f>{x)<f>,,{x) 

+  -~rc,-u+  — t - — 

C2  2c4  [<f>’(x)]2 


X 


:uuH - 


(5) 


Now,  for  lattice-Boltzmann  models  on  a  Bravais  lattice,  it  is 
well  known  that  a  Chapman-Enskog  analysis  based  on  the 
equilibrium  distribution 


where  C  and  a  are  arbitrary  constants,  and  y  is  to  be  deter¬ 
mined.  We  quickly  find  that  y  must  be  either  0  or  -D/2. 
Since  a  constant  <f>  would  not  yield  a  well  defined  h  ' ,  we  see 
that  we  must  have  <f>(x)  =  CDi\x-aCYDa,  whence  hf(x) 
=  C(a  +  x_2/D),  and  this  integrates  to  give 


AW* 


f  h0+C[ax  +  \nx ] 

h0  +  C 

r+(  1-2/D  jj 

if  D  =  2 
if  D^2, 


(7) 


where  h0  is  constant.  In  fact,  the  only  effect  of  nonzero  h0  is 
to  introduce  an  additive  constant  to  H ,  and  the  only  effect  of 
nonunity  C  is  to  scale  H  by  a  constant  factor.  In  other  words, 
h(x)  is  uniquely  specified  only  to  within  additive  and  mul¬ 
tiplicative  constants.  With  this  understanding,  we  may  say 
that  the  requirement  of  Galilean  invariance  has  uniquely 
specified  the  choice  of  H.  We  also  note  that  lim^oA'C.*:) 
=  00,  so  the  non-negativity  constraint  will  be  enforced  by  the 
dynamics. 

Finally,  we  write  the  global  Lyapunov  function  H=^XH 
by  summing  /i(Af,(x,r))  over  the  lattice.  Since  the  total  mass 
is  conserved  we  have  complete  freedom  to  choose  a,  and  so 
to  within  additive  and  multiplicative  constants  H  may  be 
written 


H(t)  oc 


2  2  ln[N,(x,f)]  for  D  =  2 


^  X'  [Ni(x,f)]‘  VD-Ni(x,t) 
x  i  2  ID 


for  D*  2, 


N'q=  — 
1  mb 


D  D(D  +  2) 

1  +  ^C'-U+  ~^~8\ C|C'“  7)  1j:uu+  •  •  • 


(6) 

will  give  rise  to  the  incompressible  Navier-Stokes  equations 
Vu=0 
and 

du  1  „  - 

—  +  gu-  Vu= - VP+  ^V2u. 

dt  6  p 


Comparing  Eqs.  (5)  and  (6),  we  identify 


for  appropriate  choices  of  a  and  C .  This  has  the  form  of  a 
Burg  entropy  [11]  for  D  =  2,  and  a  subadditive  Tsallis  en¬ 
tropy  [13]  with  parameter 

2 

for  D^2.  We  note  that  D^2  corresponds  to  q^0y  and  D 
>2  corresponds  to  0<q<  1 .  It  is  interesting  that  it  is  only  in 
the  infinite-dimensional  limit,  D— >°°,  where  the  set  of  ve¬ 
locities  becomes  infinite,  that  q~+ 1  and  we  recover  the 
Boltzmann-Gibbs  entropy  [14].  We  might  also  expect  the 
limit  of  large  b  at  constant  D  to  yield  the  Boltzmann-Gibbs 
entropy,  but  that  demonstration  will  require  the  multispeed 
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generalization  of  the  present  analysis;  that  work,  which  will 
also  provide  the  details  of  the  Chapman-Enskog  analysis,  is 
in  progress  [12].  The  numerical  implementation  of  the  model 
described  herein  is  likely  to  require  some  careful  algorithmic 
optimization,  and  is  likewise  left  to  future  publication. 

The  appearance  of  the  Burg  and  Tsallis  entropies  in  this 
context  is  fascinating.  In  a  footnote  of  their  recent  review, 
Succi,  Karlin,  and  Chen  [10]  noted  that  the  entropy  that  gave 
rise  to  the  above-mentioned  solvable  model  for  a  compress¬ 
ible  fluid  was  related  to  the  Tsallis  entropy  with  <7  =  3/2,  so 
there  may  be  more  than  one  connection  with  Tsallis  thermo¬ 
statistics  [13]  lurking  here.  There  are  precious  few  situations 
in  which  the  origins  of  Tsallis  thermostatistics  can  be  traced 
analytically  to  an  underlying  microscopic  dynamical  model, 
as  we  have  done  here. 

IV.  CONCLUSIONS 

We  have  presented  Galilean-invariant,  entropic  lattice- 
Boltzmann  models  for  the  incompressible  Navier-Stokes 
equations.  We  expect  that  these  models  will  be  useful  for  the 
simulation  of  two-  and  three-dimensional  turbulence.  As 
noted  by  Succi,  Karlin,  and  Chen  [10],  the  problem  of  find- 
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ing  perfect  models  for  lattice  models  of  the  compressible 
Navier-Stokes  equations  is  much  more  difficult  and  may  well 
be  impossible.  We  found  it  interesting  that  the  simpler  prob¬ 
lem,  for  incompressible  fluids,  is  itself  very  nontrivial.  In 
particular,  the  appearance  of  the  Burg  and  Tsallis  entropies 
for  the  H  function  is  surprising.  These  entropies  have  here¬ 
tofore  been  associated  with  long-range  interactions,  long¬ 
time  memory,  or  a  fractal  space-time  structure.  This  work 
indicates  that  they  may  also  be  relevant  to  models  with  dis¬ 
cretized  space-time  and  finite  domain  connectivity,  and  this 
surely  warrants  future  study. 
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